ZhongXu has asked that I make him the HODs from the SHAMs I made him. This should be pretty straightforward so I'll do it quickly here.

GIVING UP, moving to ds14b


In [60]:
import numpy as np
import astropy
from pearce.mocks import cat_dict
from pearce.mocks.assembias_models.table_utils import compute_prim_haloprop_bins
import h5py
from os import path

In [61]:
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
from itertools import cycle

In [62]:
config_fname = '/home/users/swmclau2/Git/pearce/bin/mcmc/nh_gg_sham_hsab_mcmc_config.yaml'

In [63]:
import yaml
with open(config_fname, 'r') as ymlfile:
    cfg = yaml.load(ymlfile)['data']

In [64]:
sim_cfg = cfg['sim']
obs_cfg = cfg['obs']
cov_cfg = cfg['cov']

In [65]:
#sim_cfg['simname'] = 'ds_14_b_sub'
#sim_cfg['halo_property'] = 'halo_vmax@mpeak'
sim_cfg['gal_property_fname'] = '/scratch/users/swmclau2/smf_dr72bright34_m7_lowm.dat'

In [66]:
cat = cat_dict[sim_cfg['simname']](**sim_cfg['sim_hps'])  # construct the specified catalog!

# TODO logspace
r_bins = obs_cfg['rbins']

obs = obs_cfg['obs']

if type(obs) is str:
    obs = [obs]

meas_cov_fname = cov_cfg['meas_cov_fname']
emu_cov_fname = cov_cfg['emu_cov_fname']
if type(emu_cov_fname) is str:
    emu_cov_fname = [emu_cov_fname]

assert len(obs) == len(emu_cov_fname), "Emu cov not same length as obs!"

n_bins = len(r_bins)-1
data = np.zeros((len(obs)*n_bins))
assert path.isfile(meas_cov_fname), "Invalid meas cov file specified"
try:
    cov = np.loadtxt(meas_cov_fname)
except ValueError:
    cov = np.load(meas_cov_fname)

assert cov.shape == (len(obs)*n_bins, len(obs)*n_bins), "Invalid meas cov shape."

In [80]:
cat.filenames[-1]


Out[80]:
'/scratch/users/swmclau2/hlists/MDHR/hlist_1.00110.list'

In [82]:
from halotools.sim_manager import RockstarHlistReader

In [86]:
reader = RockstarHlistReader(cat.filenames[-1], cat.columns_to_keep, cat.cache_filenames[-1],\
                             cat.simname, 'rockstar', 0.0,  cat.version_name, cat.Lbox, cat.pmass,\
                             overwrite=False, header_char = '#')



The information about your ascii file and the metadata about the catalog 
have been processed and no exceptions were raised. 
Use the ``read_halocat`` method to read the ascii data, 
setting the write_to_disk and update_cache_log arguments as you like. 
See the docstring of the ``read_halocat`` method
for details about these options. 


In [87]:
reader.read_halocat(cat.columns_to_convert)


...Processing ASCII data of file: 
/scratch/users/swmclau2/hlists/MDHR/hlist_1.00110.list
 
Total number of rows in detected data = 16438267
Number of rows in detected header = 48 

... working on chunk 0 of 18
... working on chunk 1 of 18
... working on chunk 2 of 18
... working on chunk 3 of 18
... working on chunk 4 of 18
... working on chunk 5 of 18
... working on chunk 6 of 18
... working on chunk 7 of 18
... working on chunk 8 of 18
... working on chunk 9 of 18
... working on chunk 10 of 18
... working on chunk 11 of 18
... working on chunk 12 of 18
... working on chunk 13 of 18
... working on chunk 14 of 18
... working on chunk 15 of 18
... working on chunk 16 of 18
... working on chunk 17 of 18
Total runtime to read in ASCII = 12.9 minutes



In [89]:
reader.halo_table


Out[89]:
<Table length=16438267>
halo_upidhalo_vacchalo_vmaxhalo_rs_klypinhalo_snapnumhalo_macchalo_halfmass_scalehalo_yhalo_idhalo_xhalo_vxhalo_vyhalo_vzhalo_m200bhalo_rshalo_rvirhalo_vpeakhalo_zhalo_mvirhalo_nfw_conchalo_hostid
int64float32float32float32int64float32float32float32int64float32float32float32float32float32float32float32float32float32float32float32int64
-1210300.01522.560.48661742159500.0108100.0478.52773330340320.6493282.17815.03277.491.169e+150.5128372.09695137900.0688.6581.041e+154.08893733303403
-190930.01447.710.54400842182000.0115600.0328.863733018090119.928-621.74-772.24333.821.0083e+150.5931022.02309112800.0612.1999.351e+143.41103733018090
-1142800.01365.170.73713542345300.099720.0290.43573303497034.4806111.96-487.09-27.311.0533e+150.7529881.9521163760.0660.0228.401e+142.59249733034970
-146020.01464.010.2187954249260.041620.0389.7277333035510.43736-130.27-368.0268.377.7247e+140.2288171.814647730.0655.7486.748e+147.93037733303551
-1-2483.01360.740.27949342-23930.074800.0485.214735444141171.419-274.79-200.48478.967.2269e+140.3097541.7730287440.0688.5376.295e+145.72395735444141
-116020.01294.130.3405854215480.048080.0330.24873303515194.7808-275.94-170.96-325.046.926e+140.346831.7428344370.0641.9355.979e+145.02503733035151
-142640.01275.640.37478742-9852.062450.0290.60773301794662.462630.09-1.4168.156.5839e+140.3681251.7411285280.0623.4125.961e+144.72969733017946
-1127900.01265.650.3402234218770.060780.0253.68373303523044.2488196.94233.4545.586.385e+140.3880691.7095781810.0654.6225.643e+144.40532733035230
-181700.01291.390.2850354224870.063390.0468.582735444356184.141-398.99199.15464.135.6973e+140.3243841.6995688690.0685.3085.544e+145.23935735444356
-139570.01237.160.3776424223930.039080.0491.44735427805134.029-114.0227.2198.026.2301e+140.3616751.6957845940.0591.8475.507e+144.6887735427805
...............................................................
-10.038.580.0135823420.00.0747.327749194485911.659315.3-106.03-214.783.4892e+100.0135820.0536560.0900.5391.745e+103.95052749194485
-1-4.34438.570.0136374420.0-1.448747.223749206399946.27543.86-5.89-96.686.1061e+100.0136370.053656-2.172939.5391.745e+103.93459749206399
-1-8.68741.480.0081129242-40.91-2.896748.624749208090950.429-32.37-78.68-127.639.5953e+100.0081130.053656-4.344933.4151.745e+106.61358749208090
-10.050.820.003285420.00.0748.072749208624955.877180.13-249.8-42.811.7446e+100.0032850.0536560.0972.4431.745e+1016.3336749208624
-1-4.34437.380.0244003420.0-5.793748.92749209830960.318154.2641.8953.927.8507e+100.02440.053656-8.69992.7191.745e+102.19902749209830
-1-13.8643.850.00607511420.05.449748.866749210253954.215143.62-125.2732.891.3085e+110.0060750.053656-5.609988.251.745e+108.83226749210253
-10.044.570.00562847420.00.0748.487749210436899.744100.92-227.97-65.421.134e+110.0056280.0536560.0954.7121.745e+109.53376749210436
-10.039.630.0108573420.00.0747.821749210819955.974179.3-226.36-80.179.5953e+100.0108570.0536560.0973.9021.745e+104.94206749210819
-1-4.34439.630.010852842-40.91-1.448747.647749210821957.03490.79-202.35-63.889.5953e+100.0108530.053656-2.172968.6751.745e+104.94389749210821
-10.038.060.0158003420.00.0746.808749211177981.854182.63-45.34-44.33.4892e+100.01580.0536560.0928.5311.745e+103.39595749211177

In [67]:
#cat.load(sim_cfg['scale_factor'], **sim_cfg['sim_hps'])
#cat.populate()# will generate a mock for us to overwrite
gal_property = np.loadtxt(sim_cfg['gal_property_fname'])
halo_property_name = sim_cfg['halo_property']
min_ptcl = sim_cfg.get('min_ptcl', 200)
nd = float(sim_cfg['nd'])
scatter = float(sim_cfg['scatter'])

In [91]:
cat.h


Out[91]:
0.7020000000000001

In [90]:
from AbundanceMatching import *
af =  AbundanceFunction(gal_property[:,0], gal_property[:,1], sim_cfg['af_hyps'], faint_end_first = sim_cfg['reverse'])
remainder = af.deconvolute(scatter, 20)
# apply min mass
halo_table = reader.halo_table#cat.halocat.halo_table#[cat.halocat.halo_table['halo_mvir']>min_ptcl*cat.pmass] 
nd_halos = calc_number_densities(halo_table[halo_property_name], cat.Lbox) #don't think this matters which one i choose here
catalog_w_nan = af.match(nd_halos, scatter)
n_obj_needed = int(nd*(cat.Lbox**3))
catalog = halo_table[~np.isnan(catalog_w_nan)]
sort_idxs = np.argsort(catalog)

In [69]:
np.all(cat.halocat.halo_table['halo_upid']==-1)


Out[69]:
False

In [92]:
np.sum(halo_table['halo_upid'] == -1)


Out[92]:
14565913

In [93]:
final_catalog = halo_table[~np.isnan(catalog_w_nan)][sort_idxs[:n_obj_needed]]

#final_catalog['x'] = final_catalog['halo_x']
#final_catalog['y'] = final_catalog['halo_y']
#final_catalog['z'] = final_catalog['halo_z']
#final_catalog['halo_upid'] = -1
# FYI cursed.
#cat.model.mock.galaxy_table = final_catalog
# TODO save sham hod "truth"
for idx, (o, ecf) in enumerate(zip(obs, emu_cov_fname)): calc_observable = getattr(cat, 'calc_%s' % o) y = calc_observable(r_bins) data[idx*n_bins:(idx+1)*n_bins] = y

In [94]:
final_catalog.colnames


Out[94]:
['halo_upid',
 'halo_vacc',
 'halo_vmax',
 'halo_rs_klypin',
 'halo_snapnum',
 'halo_macc',
 'halo_halfmass_scale',
 'halo_y',
 'halo_id',
 'halo_x',
 'halo_vx',
 'halo_vy',
 'halo_vz',
 'halo_m200b',
 'halo_rs',
 'halo_rvir',
 'halo_vpeak',
 'halo_z',
 'halo_mvir',
 'halo_nfw_conc',
 'halo_hostid']

In [95]:
from halotools.mock_observables import *

In [96]:
x, y, z = [final_catalog[c] for c in ['halo_x', 'halo_y', 'halo_z']]
pos = return_xyz_formatted_array(x, y, z, period=cat.Lbox)

In [97]:
xi_all = tpcf(pos / cat.h, r_bins, period=cat.Lbox / cat.h, num_threads=4,
                              estimator='Landy-Szalay')

In [98]:
xi_all


Out[98]:
array([ -1.00000000e+00,  -1.00000000e+00,   2.34036457e+01,
         1.52151441e+02,   2.97679025e+02,   2.91343348e+02,
         1.93159382e+02,   1.09047017e+02,   5.26977826e+01,
         2.62023209e+01,   1.34719705e+01,   7.02559721e+00,
         3.58227506e+00,   1.87404796e+00,   1.03531472e+00,
         5.79425318e-01,   3.12589004e-01,   1.57569414e-01])

In [100]:
plt.plot(xi_all)
plt.yscale('log')



In [102]:
import h5py
f = h5py.File('/home/users/swmclau2/scratch/PearceMCMC/pearce_mcmc_nh_gg_sham_hsab.hdf5', 'r')

In [103]:
cov = f['cov'].value

In [105]:
np.savetxt('/home/users/swmclau2/Git/pearce/bin/shams/xigg_cov_mcmc.npy', cov)

In [ ]: